surge analysis on IOC tide gauge¶

case study: cyclone FUNG WONG

Install¶

Libraries needed for this notebook:

pip install searvey hvplot utide ipykernel geoviews pyarrow

In [31]:
import searvey
import utide
import pandas as pd
import hvplot.pandas

from shapely.geometry import box
In [33]:
ioc_stations = searvey.get_ioc_stations()
ioc_stations
Out[33]:
ioc_code gloss_id lat lon country location connection contacts dcp_id last_observation_level ... observations_expected_per_week observations_ratio_per_week observations_arrived_per_month observations_expected_per_month observations_ratio_per_month observations_ratio_per_day sample_interval average_delay_per_day transmit_interval geometry
0 aarh <NA> 56.150 10.220 Denmark Aarhus ftp Danish Meteorological Institute ( Denmark ) NaN 0.26 ... 1008.0 99 4407 4464.0 99 97 10' 6' 10' POINT (10.22 56.15)
1 abas 327 44.020 144.290 Japan Abashiri SWJP40 Japan Meteorological Agency ( Japan ) ABASHIRI 2.10 ... 10080.0 100 43970 44640.0 98 98 1' 8' 10' POINT (144.29 44.02)
2 abed <NA> 57.140 -2.080 UK Aberdeen ftp National Oceanography Centre ( UK ) NaN 3.25 ... 672.0 100 2975 2976.0 100 100 15' 28' 15' POINT (-2.08 57.14)
3 abur 82 31.580 131.410 Japan Aburatsu SWJP40 Japan Meteorological Agency ( Japan ) ABURATSU 2.50 ... 10080.0 100 43970 44640.0 98 98 1' 8' 10' POINT (131.41 31.58)
4 acaj 182 13.574 -89.838 El Salvador Acajutla SV SZXX01 Ministerio de Medio Ambiente y Recursos Natura... 300434064008810 1.16 ... 10080.0 91 36760 44640.0 82 90 1' 6' 5' POINT (-89.838 13.574)
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1720 zkth <NA> 37.781 20.905 Greece Zakynthos, Ionian Islands bgan National Observatory of Athens ( Greece ) GR-ZKTH-00 0.01 ... 10080.0 0 NaN NaN 0 0 1' NaN 1' POINT (20.905 37.781)
1721 zldw <NA> 51.346 3.200 Belgium Zeebrugge Leopold II dam web Maritieme Dienstverlening en Kust ( Belgium ) NaN 3.86 ... 2016.0 100 8902 8928.0 100 100 5' 8' 5' POINT (3.2 51.346)
1722 zwdw <NA> 51.355 3.178 Belgium Zeebrugge Wielingendok web Maritieme Dienstverlening en Kust ( Belgium ) NaN 3.88 ... 2016.0 100 8908 8928.0 100 100 5' 8' 5' POINT (3.178 51.355)
1723 zygi <NA> 34.727 33.338 Cyprus Zygi ftp Cyprus Oceanography Center ( Cyprus ) NaN 1.97 ... 20160.0 0 -down- NaN 0 0 0.5' NaN 1' POINT (33.338 34.727)
1724 zygi1 <NA> 34.726 33.340 Cyprus Zygi bgan Cyprus Marine and Maritime Institute ( Cyprus ... ZYGI1 0.63 ... 10080.0 0 -down- NaN 0 0 1' NaN 1' POINT (33.34 34.726)

1725 rows × 25 columns

In [34]:
bbox = box(100.0, -5.0, 130.0, 25.0)  # lon_min, lat_min, lon_max, lat_max
In [36]:
selected_stations = ioc_stations[ioc_stations.within(bbox)]
selected_stations
Out[36]:
ioc_code gloss_id lat lon country location connection contacts dcp_id last_observation_level ... observations_expected_per_week observations_ratio_per_week observations_arrived_per_month observations_expected_per_month observations_ratio_per_month observations_ratio_per_day sample_interval average_delay_per_day transmit_interval geometry
45 ambon 68 -3.683 128.183 Indonesia Ambon SZXX01 Geospacial Agency of Indonesia ( Indonesia ) +... 301434060409970 4.46 ... 10080.0 91 36510 44640.0 82 90 1' 5' 5' POINT (128.183 -3.683)
47 AMTSI <NA> -0.847 121.601 Indonesia Port of Ampana, Central Sulawesi web Kepala Badan Meteorologi, Klimatologi, dan Geo... NaN 0.01 ... 10080.0 96 37857 44640.0 85 98 1' 2' 5' POINT (121.601 -0.847)
132 BETSI <NA> 4.230 126.788 Indonesia Port of Beo, North Sulawesi web Kepala Badan Meteorologi, Klimatologi, dan Geo... NaN 0.01 ... 10080.0 96 43470 44640.0 97 76 1' 1' 5' POINT (126.788 4.23)
136 bitu 69 1.439 125.190 Indonesia Bitung SZXX01 Geospacial Agency of Indonesia ( Indonesia ) +... 300434067058990 8.15 ... 10080.0 90 36675 44640.0 82 89 1' 5' 5' POINT (125.19 1.439)
183 bupo <NA> -1.030 100.396 Indonesia Bungus Port bgan Joint Research Centre ( Europe ) +Ministry of ... BUPO 1.45 ... 10080.0 0 -down- NaN 0 0 1' NaN 1' POINT (100.396 -1.03)
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1621 txil <NA> 22.353 120.383 Taiwan Xiao Liuqiu web Central Weather Administration ( Taiwan ) NaN 0.20 ... 168.0 100 743 744.0 100 100 60' 16' 1h POINT (120.383 22.353)
1622 tyon <NA> 22.819 120.198 Taiwan Yongan web Central Weather Administration ( Taiwan ) NaN -999.00 ... 168.0 100 743 744.0 100 100 60' 16' 1h POINT (120.198 22.819)
1672 vung <NA> 10.340 107.071 Viet Nam Vung Tau SZXX01 Hydro-meteorological and Environmental Station... 300434067059990 2.26 ... 10080.0 88 35835 44640.0 80 87 1' 5' 6' POINT (107.071 10.34)
1706 xish <NA> 16.830 112.330 China Xisha SZCI01 China Meteorological Administration ( China ) 11745 1.30 ... 10080.0 0 -down- NaN 0 0 1' NaN 5' POINT (112.33 16.83)
1717 zhap 78 21.580 111.820 China Zhapo SZCI01 China Meteorological Administration ( China ) 09731 2.42 ... 10080.0 44 18491 44640.0 41 45 1' 7' 1' POINT (111.82 21.58)

70 rows × 25 columns

In [32]:
plot_ = selected_stations.hvplot.points(
    geo=True,
    tiles=True,
    hover_cols=["ioc_code"],
    s=100,
    line_color='k',
    title = "stations in the selected region"
).opts(width=600, height=700)
plot_
Out[32]:

detide the stations¶

In [30]:
def surge(ts: pd.Series, lat: float, rsmp: int = None):
    ts0 = ts.copy()
    OPTS = {
        "constit": "auto",
        "method": "ols", # ols is faster and good for missing data (Ponchaut et al., 2001)
        "order_constit": "frequency",
        "Rayleigh_min": 0.97,
        "lat": lat,
        "verbose": True,
    }
    if rsmp is not None:
        ts = ts.resample(f"{rsmp}min").mean()
        ts = ts.shift(freq=f"{rsmp / 2}min")
    coef = utide.solve(ts.index, ts, **OPTS)
    tidal = utide.reconstruct(ts0.index, coef, verbose = OPTS['verbose'])
    return pd.Series(data=ts0.values - tidal.h, index=ts0.index)

download the data¶

create a data folder with raw and surge subfolders

In [12]:
! mkdir -p data/{raw,surge}

let's first download all the station in the bounding box

In [16]:
drop_columns = ["bat"]

def serialize(d): # to export metadata
    out = {}
    for k, v in d.items():
        if v is not pd.NA and not pd.isna(v):
            out[k] = str(v)
    return out
In [17]:
for irow, row in selected_stations.iterrows():
    lat = row.lat
    ioc_code = row.ioc_code
    if ioc_code in list(detide_selection.keys()):
        df_raw = searvey.fetch_ioc_station(
            ioc_code,
            pd.Timestamp.now()-pd.Timedelta(days=90), # we need at least 90 days (in theory we'd need more..) to remove properly tidal constituents
            pd.Timestamp.now()
        )
        df_raw.attrs = {**serialize(dict(row)), **{"signal_type": "raw"}}
        df_raw.to_parquet(f"data/raw/{ioc_code}.parquet")
In [6]:
for irow, row in selected_stations.iterrows():
    df_raw = pd.read_parquet(f"data/raw/{row.ioc_code}.parquet")
    df_raw.loc["2025-11-01":].dropna().hvplot(
        title = f"detided signal for {row.ioc_code}, {row.location}, lat:{row.lat}, lon:{row.lon}"
    )
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
Out[6]:
bokeh backend could not plot any Elements in the Overlay.
Out[6]:
:NdOverlay   [Variable]
Out[6]:

detide the station¶

Let's look at the stations of interest

In [ ]:
In [18]:
detide_selection = {
    "lega" :"rad",
    "luba" :"prs",
    "quin" :"ras",
    "qing" :"rad",
    "quar" :"flt",
    "mani" :"prs",
    "subi" :"rad",
    "thsi" :"rad",
    "txil" :"rad",}
In [27]:
for irow, row in selected_stations.iterrows():
    if row.ioc_code in list(detide_selection.keys()):
        df_raw = pd.read_parquet(f"data/raw/{row.ioc_code}.parquet")
        df_surge = surge(df_raw[detide_selection[row.ioc_code]].dropna(), lat=row.lat, rsmp=2)
        df_surge.loc["2025-11-01":].hvplot(
            title = f"detided signal for {row.ioc_code}, {row.location}, lat:{row.lat}, lon:{row.lon}"
        )
        df_surge.attrs = dict(row)
        df_surge.attrs = {**serialize(dict(row)), **{"signal_type": "detide"}}
        df_surge.to_frame().to_parquet(f"data/surge/{row.ioc_code}.parquet")

solve: matrix prep ... solution ... done.
prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done.
prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done.
prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done.
prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done.
prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done.
prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done.
prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done.
prep/calcs ... done.
Out[27]:
solve: matrix prep ... solution ... done.
prep/calcs ... done.
Out[27]:
In [26]:
(plot_*ioc_stations[ioc_stations.ioc_code.isin(list(detide_selection.keys()))].hvplot.points(
    geo=True,
    tiles=True,
    c = 'r',
    s=50,
    hover_cols=["ioc_code"],
    title = "stations that recorded a surge"
)).opts(width=800, height=800)
Out[26]:

look at a specific station¶

In [ ]:
station = "lega"
for irow, row in selected_stations.iterrows():
    if row.ioc_code == station:
        df_raw = pd.read_parquet(f"data/raw/{row.ioc_code}.parquet")
        df_raw
        station
        df_raw.loc["2025--01":].dropna().hvplot(
            title = f"signal for {row.ioc_code}, {row.location}, lat:{row.lat}, lon:{row.lon}"
        )